Joint estimation of survival and dispersal effectively corrects the permanent emigration bias in mark-recapture analyses

Robust and reliable estimates of demographic parameters are essential to understand population dynamics. Natal dispersal is a common process in monitored populations and can cause underestimations of survival and dispersal due to permanent emigration. Here, we present a multistate Bayesian capture-mark-recapture approach based on a joint estimation of natal dispersal kernel and detection probabilities to address biases in survival, dispersal, and related demographic parameters when dispersal information is limited. We implement this approach to long-term data of a threatened population: the Bonelli’s eagle in Catalonia (SW Europe). To assess the method’s performance, we compare demographic estimates structured by sex, age, and breeding status in cases of limited versus large data scales, with those of classical models where dispersal and detection probabilities are estimated separately. Results show substantial corrections of demographic estimates. Natal dispersal and permanent emigration probabilities were larger in females, and consequently, female non-breeder survival showed larger differences between separate and joint estimation models. Moreover, our results suggest that estimates are sensitive to the choice of the dispersal kernel, fat-tailed kernels providing larger values in cases of data limitation. This study provides a general multistate framework to model demographic parameters while correcting permanent emigration biases caused by natal dispersal.

. State transition matrix of the multistate joint-estimation models. Rows indicate state of individual i at time t, columns indicate state of individual i at time t+1. The matrix includes the parameters φNB (nonbreeder survival), φB (breeder survival), (recruitment to the breeding population), Pd (joint-estimation parameter), and r (probability of encountering a dead marked individual that was either non-breeder, rNB, or breeder, rB). Subindices indicate age (a), sex (s), individual (i), and time (t). States are Alive as Non-Breeder (ANB), Alive as Breeder (AB), Dead as Non-Breeder (DNB), Dead as Breeder (DB), and Absorbing State (AS). Due to the length of the first absorbing state expression (i.e., transition from ANB in t to AS in t+1), we wrote this expression outside the matrix (*A). Note that r (recovery) is indeed an observation parameter and is often modelled in the observation matrix, but for convergence issues in Bayesian models it must be modelled in the state matrix. Refer to Kéry and Schaub (2012) for further information.  Figure S2. Observation matrix of the multistate models. Rows indicate state of individual i at time t, columns indicate observation/event of individual i at time t. The matrix includes the parameters PNB (non-breeder recapture), PB1 (probability of being observed for the first time once the individual has become breeder), and PB (breeder recapture for those individuals observed as breeders at least once already), and the binary covariate covpb, which assigns either PB1 or PB as observation parameters to deal with detection heterogeneity in breeder detection probabilities.     Table S2. Long-distance natal dispersal probability (i.e., distances higher than 200 km) by dispersal kernel and sex.

Supplementary Material S3: A brief guide on adapting the Joint Estimation multistate matrices to different study systems
The joint estimation approach is a flexible formulation to address biases in demographic parameters caused by permanent emigration. This model formulation is divided into a spatial submodel and a multistate capture-mark-recapture (CMR) submodel. The spatial submodel models natal dispersal, while the CMR submodel provides information on demographic parameters such as survival and sexual maturity (i.e., age of first breeding or recruitment to the breeding population, among other possibilities). Both submodels are linked by a specific parameter that estimates individual probabilities of permanent emigration accounting for information provided by both submodels and detectability along space. This parameter is accommodated in the state matrices of the multistate submodel to represent the permanent emigration process. Whenever an individual achieves sexual maturity, it undertakes natal dispersal towards the area of breeding. If that site is outside the study area, this is considered as permanent emigration. Such process is modelled as a transition to the absorbing state (i.e., a state of permanent undetectability) as the individual at issue will not be observable anymore throughout the course of the study. Based on this modelling approach, to implement the joint estimation it is only necessary to model survival and sexual maturity in the state matrix of the multistate submodel ( Figure S6). As for the observation matrix, a single recapture probability would be necessary, same for breeders and non-breeders ( Figure S7). Note that except for the Pd parameter, which should be structured by individual in case we want to account for birth territory locations and the border effect, structuring parameters by age, sex, or sexual maturity status as in our study, is not necessary either for implementing this method.
The flexibility of the multistate capture-mark-recapture method may facilitate the adaptation of the joint estimation approach to other study designs. For instance, note that in Figure S6 we modelled the non-breeder to breeder transition as a stochastic probability, but if individuals become sexually mature at a fixed age (i.e., deterministically), this may also be adapted in the matrices using age states. As an example of an appropriate formulation, if an individual becomes sexually mature at the end of its second year of life, one may use the states "Alive -Non Breeder 1 year old", "Alive -Non Breeder 2 years old", "Alive -Breeder", and the Absorbing State (Figures S8 and S9). In this case, transitions from live states should only be determined by survival, as the probability to transition from being non-breeder at year 1 to being non-breeder at year 2, and to being breeder at year 3, only depends on surviving from year to year.
From these basic models, it is easy to add complexity to the state transition and observation matrices in further steps. For instance, from the simplest matrices presented here ( Figures S6   and S7), we may add the modelling of recoveries of dead individuals (Figures S10 and S11.
Recall that while recoveries have traditionally been modelled as observation parameters, in Bayesian hierarchical models and BUGS language we usually model them in the state matrices to avoid convergence issues.

ANB-2Y
As shown, using the flexibility of multistate matrix we can add further layers of complexity to basic matrices in order to accommodate study systems with greater structural or modelling particularities. In the particular case of our study, we easily adapted the matrices to the particularities of the study species, with detection heterogeneity in breeding birds, and recovery probabilities. Further specificities of study populations may be added, such as disease status [1,2], body condition [3], and tag loss [4], among others. The prior distributions of the mean and standard deviation in the lognormal distribution and shape and rate in the gamma distribution will determine the shape of the prior dispersal kernels.
Recall that in our study both distributions are truncated at a maximum value of 1200 km. A fundamental step is to find prior distributions that provide similar dispersal distance expectations for both models. This way, we can ensure priors do not affect the estimates of our models based on different distributions in significantly different ways. The information provided by prior beliefs in a model can be easily evaluated using prior predictive checks (PPC; [2,3]). PPCs are based on simulating and visualizing data from prior distributions. Using PPCs, our aim was to choose weakly informative priors that 1) encompass the range of potential values of each parameter, 2) generate reasonable values of natal dispersal distances, that is, those that encompass the range of observed distances for the Bonelli's eagle in continental western Europe in published studies [4], and 3) regularize against dispersal distance values that are considered infrequent or have never been reported according to the literature [2,4,5].    [4].
When compared to the published data in the study population and the neighbouring French population, the means of our prior dispersal kernels for both distributions are lower ( Figure   S13). This is because 1) dispersal distances on published articles were analysed using normal distributions, which have markedly different shapes to gamma and lognormal distributions, and 2) while published data do not account for sex, we will be splitting our dispersal estimations into sexes, and males presumably have shorter dispersal distances than females, thus we prefer to be conservative. One further possibility may be to use priors that reflect wider uncertainties about the dispersal kernels of both distributions between the range of distances considered in our study (i.e., 0-1200 km). Although this may also be an interesting choice, prior dispersal kernel distributions would then keep considerable prior probability densities at large and unlikely distances, which may not contribute to regularizing towards more likely values and thus not be as optimal as the priors selected in this study. However, an interesting exercise to understand the effect of our prior selection would be to analyse the sensitivity of prior choice on posterior inference (i.e., narrower weakly regularizing prior as the one chosen in our study and shown in Figure S12 vs. wider prior shown in Figure S14). We selected such wider priors by trial and error and using PPCs. We used normally-distributed priors truncated at 0 for each distribution, again focusing on obtaining similar prior distance expectations (lognormal and gamma respectively). These alternative dispersal kernels show median distances of 212 and 202 km, 2.5% quantile distances of 17 and 16 km, and 97.5% quantile distances of 1080 and 1076 km for the lognormal and the gamma distribution respectivelythat is, both kernel distributions are very similar, and show way larger expectations than the ones selected in this study ( Figure S14).
We compared posterior estimates of our models when using both sets of priors [1]. We did so for SEP-CAT and JOINT-CAT models (we excluded JOINT-ALL models because of their longer computation times). We compared the survival and average dispersal distance estimates obtained using both priors in the range of different models considered in this sensitivity analysis. We define the prior selected in our study as "Narrower" and the alternative prior described above as "Wider". The sensitivity analysis reveals very mild differences in the posterior estimates of survival ( Figure S15) and average dispersal distance ( Figure S16) when using the narrower priors selected in this study and alternative wider priors. Estimates of both survival and dispersal generally show very limited increasesif any -when using wider prior dispersal kernels with larger expectancies at large distances. Such slight increases appear to occur in both distributions. This suggests that 1) the effect of our selected priors on posterior inference is limited, and 2) the differences found between models and distributions are not explained by our prior dispersal kernel choices. Figure S16. Posterior estimates of average dispersal distance (km) of SEP-CAT and JOINT-CAT models when using narrower priors (i.e., the ones selected in this study) with alternative wider priors. Circles and squares indicate median values. Thick lines indicate the 66% HPDI and fine lines show the 89% HPDI.